In [1]:

    

from pearce.mocks import cat_dict
import numpy as np
from os import path


    

In [2]:

    

import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()


    

In [3]:

    

a = 1.0#0.81120
z = 1.0/a - 1.0


    

In [4]:

    

print z


    

    




0.0

In [5]:

    

cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[a]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!

cat.load_catalog(a)
#halo_masses = cat.halocat.halo_table['halo_mvir']


    

In [6]:

    

cat.load_model(a, 'redMagic')


    

In [7]:

    

params = cat.model.param_dict.copy()
params['mean_occupation_centrals_assembias_param1'] = 0.0
params['mean_occupation_satellites_assembias_param1'] = 0.0
params['logMmin'] = 12.089
params['sigma_logM'] = 0.33
params['f_c'] = 1.0
params['alpha'] = 1.1
params['logM1'] = 13.3
params['logM0'] = params['logMmin']

print params


    

    




{'logM1': 13.3, 'mean_occupation_satellites_assembias_param1': 0.0, 'logMmin': 12.089, 'mean_occupation_centrals_assembias_param1': 0.0, 'f_c': 1.0, 'logM0': 12.089, 'sigma_logM': 0.33, 'alpha': 1.1}

In [8]:

    

cat.populate(params)


    

In [9]:

    

from halotools.empirical_models import PrebuiltSubhaloModelFactory
model = PrebuiltSubhaloModelFactory('behroozi10')

from halotools.sim_manager import CachedHaloCatalog
halocat = CachedHaloCatalog(simname = 'bolshoi', redshift=0, version_name ='halotools_alpha_version2')
model.populate_mock(halocat)


    

In [10]:

    

#theta_bins = np.logspace(np.log10(0.004), 0, 24)#/60
#tpoints = (theta_bins[1:]+theta_bins[:-1])/2
theta_bins = np.logspace(-2,0,15)
tpoints = (theta_bins[:-1]+theta_bins[1:])/2.0


    

In [ ]:

    

from halotools.mock_observables import mock_survey, angular_tpcf


    

In [ ]:

    

#This is what is currently in Pearce
n_cores = 'max'


pos = np.vstack([model.mock.galaxy_table[c] for c in ['x', 'y', 'z']]).T
vels = np.vstack([model.mock.galaxy_table[c] for c in ['vx', 'vy', 'vz']]).T

# TODO is the model cosmo same as the one attached to the cat?
ra, dec, _ = mock_survey.ra_dec_z(pos * model.mock.cosmology.h, vels , cosmo=model.mock.cosmology)
ang_pos = np.vstack((np.degrees(ra), np.degrees(dec))).T

n_rands = 5
rand_pos = np.random.random((pos.shape[0] * n_rands, 3)) * model.mock.Lbox#*self.h
rand_vels = np.zeros((pos.shape[0] * n_rands, 3))

rand_ra, rand_dec, rand_z = mock_survey.ra_dec_z(rand_pos * model.mock.cosmology.h, rand_vels , cosmo=model.mock.cosmology)
rand_ang_pos = np.vstack((np.degrees(rand_ra), np.degrees(rand_dec))).T

# NOTE I can transform coordinates and not have to use randoms at all. Consider?
wt_all = angular_tpcf(ang_pos, theta_bins, randoms=rand_ang_pos, num_threads=n_cores)
print wt_all


    

    




/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/halotools-0.6.dev4681-py2.7-linux-x86_64.egg/halotools/mock_observables/two_point_clustering/clustering_helpers.py:134: UserWarning: 
 `sample1` exceeds `max_sample_size` 
downsampling `sample1`...
  warn(msg)

In [ ]:

    

#now, from halotools docs 
#mask = model.mock.galaxy_table['stellar_mass'] > 10**10.5
gals = model.mock.galaxy_table#[mask]
coords = np.vstack([gals['x'], gals['y'], gals['z']]).T
vels = np.vstack([gals['vx'], gals['vy'], gals['vz']]).T

ra, dec, z = mock_survey.ra_dec_z(coords*model.mock.cosmology.h, vels, cosmo=model.mock.cosmology)
ra = np.degrees(ra)
dec = np.degrees(dec)

Nran=10**5
ran_coords = np.random.random((Nran,3))*cat.model.mock.Lbox
ran_vels = np.zeros((Nran,3))

ran_ra, ran_dec, ran_z = mock_survey.ra_dec_z(ran_coords*model.mock.cosmology.h, ran_vels, cosmo=model.mock.cosmology)
ran_ra = np.degrees(ran_ra)
ran_dec = np.degrees(ran_dec)

angular_coords = np.vstack((ra,dec)).T
ran_angular_coords = np.vstack((ran_ra,ran_dec)).T


w_theta_with_randoms = angular_tpcf(angular_coords, theta_bins, randoms=ran_angular_coords, num_threads='max')
print w_theta_with_randoms


    

In [ ]:

    

print rand_pos.shape, ran_coords.shape, pos.shape


    

In [ ]:

    

plt.plot(tpoints,1+wt_all)
plt.plot(tpoints, 1+w_theta_with_randoms)
#plt.xscale('log')
plt.loglog()
plt.xlim([1e-2, 1.0])
#plt.ylim([1e-4, 2.0])


    

In [ ]:

    

coords = np.vstack((gals['x'],gals['y'],gals['z'])).T - cat.model.mock.Lbox/2.0
vels = np.vstack((gals['vx'],gals['vy'],gals['vz'])).T

ra_init, dec_init, z = mock_survey.ra_dec_z(coords*model.mock.cosmology.h, vels, cosmo=model.mock.cosmology)

#keep a complete spherical volume
r = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
keep = r<=cat.model.mock.Lbox[0]

ra = np.degrees(ra_init[keep])
dec = np.degrees(dec_init[keep])
angular_coords = np.vstack((ra,dec)).T
w_theta = angular_tpcf(angular_coords, theta_bins, num_threads='max')

print w_theta


    

In [ ]:

    

plt.plot(tpoints,1+wt_all)
plt.plot(tpoints, 1+w_theta_with_randoms)
plt.plot(tpoints, 1+w_theta)
#plt.xscale('log')
plt.loglog()
plt.xlim([1e-2, 1.0])
#plt.ylim([1e-4, 2.0])


    

In [ ]:

    

w_theta/wt_all


    

In [ ]:

    

Nran=10**6
ran_coords = np.random.random((Nran,3))*model.mock.Lbox - model.mock.Lbox/2.0
ran_vels = np.zeros((Nran,3))

ran_ra, ran_dec, ran_z = mock_survey.ra_dec_z(ran_coords, ran_vels, cosmo=model.mock.cosmology)

#keep a complete spherical volume
r = np.sqrt(ran_coords[:,0]**2 + ran_coords[:,1]**2 + ran_coords[:,2]**2)
keep = r<model.mock.Lbox[0]/2.0

ran_ra = np.degrees(ran_ra[keep])
ran_dec = np.degrees(ran_dec[keep])
ran_angular_coords = np.vstack((ran_ra,ran_dec)).T

w_theta_with_randoms = angular_tpcf(angular_coords, theta_bins, randoms=ran_angular_coords,\
                                    num_threads='max')


    

In [ ]:

    

plt.plot(tpoints,1+wt_all)
plt.plot(tpoints, 1+w_theta_with_randoms)
plt.plot(tpoints, 1+w_theta)
#plt.xscale('log')
plt.loglog()
plt.xlim([1e-2, 1.0])
#plt.ylim([1e-4, 2.0])


    

In [ ]:

    

plt.plot(tpoints, w_theta/wt_all)
plt.xscale('log')


    

In [ ]: